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One hypothesized explanation for water's anomalies imagines the existence of a 
liquid-liquid (LL) phase transition line separating two liquid phases and termi- 
nating at a LL critical point. We simulate the classic ST2 model of water for 
times up to 1000 ns and system size up to = 729. We find that for state points 
near the LL transition line, the entire system fiips rapidly between liquid states 
of high and low density. Our finite-size scaling analysis accurately locates both 
the LL transition line and its associated LL critical point. We test the stability 
of the two liquids with respect to the crystal and find that of the 350 systems 
simulated, only 3 of them crystallize and these 3 for the relatively small system 
size N=343 while for all other simulations the incipient crystallites vanish on a 
time scales smaller than ^ 100ns. 

We perform extensive molecular dynamics (MD) simulations of ST2-water in the constant- 
temperature, constant-pressure ensemble. We equilibrate the system for 1000ns for 127 
state points in the supercooled liquid region of water. Pressure P ranges from 190 MPa to 
240 MPa, while temperature T is as low as T =230 K at high P, and 244 K at low P. We 
make 624 different simulations, 341 as long as 1000 ns, and for four system sizes, = 216 
(80 state points), 343 (75 state points), 512 (44 state points), and 729 molecules (46 state 
points). For the majority of state points studied we average our results over several (< 11) 
independent runs. We interpolate our data along isobars using the histogram reweighting 
method [l|. For P > 200 MPa, we find that the density p decreases sharply within a 
narrow temperature range, while at lower P it falls off with T continuously. This behavior 
is consistent with a discontinuous phase transition at high-P between a high-density liquid 
(HDL) and a low-density liquid (LDL) ending in a liquid-liquid (LL) critical point at lower 

p{Fig.m. 

This LL critical point was hypothesized 12| based on studies of the ST2 model, and 

nn n 

subsequently studied in detail by many others using, in addition to ST2 P, TIP5P [51, 
TIP4P 3, TIP4P-EW Q and TIP4P/2005 ^ as well as coarse-grained models 9-ll|. The 
existence of the LL critical point allows one to understand X-ray spectroscopy results 
[l^ . and explains the increasing correlation length in bulk water upon cooling as found 
experimentally 15| and the hysteresis effects 16|. Holten et al. 17|, ll8| reviews available 
experimental information and shows that the assumption of a LL critical point in supercooled 
water provides an accurate account on the experimental thermodynamic properties. 
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Abrupt changes in the global density p are related to the appearance of different local 



structures. Among various parameters describing the local structures we identify [19| and 
^3, defined in the Methods Section, as good quantities to distinguish the LDL and the HDL 
phase and the best quantities to distinguish them from ice. The average values of ip^ of the 
two phases differ by about 50%, the LDL phase being characterized by greater order in the 
second shell than in the HDL phase. 



Liu et al. j4|, using histogram reweighed Monte Carlo methods in the grand canonical 
ensemble for only one but quite large system size, find an order parameter distribution 
function consistent with a critical point belonging to the universality class of a 3 dimensional 



Ising model. In Ref. [3l| Limmer and Chandler question this result using the umbrella 
sampling method to evaluate the free energy landscape of the ST2 model near a single 
state point and for a single system size (A^ = 216). They find two minima in the free 
energy landscape: one for liquids and one for crystalline structure. They do not find a 
third minimum corresponding to the LDL and conclude that the LDL does not exist as a 
metastable state, but only as a transitional state from HDL to crystal. However, Sciortino 



et al. in 



33[ show, with an implementation for 200 > A^ > 327 of the umbrella sampling 



that guarantees very high resolution in the exploration of the free energy landscape, the 

;o the LDL state metastable with respect to the 



presence of the minimum corresponding 
crystal, reconfirming the results of Ref. J] and at variance with Ref. 3l|- To contribute 
to the discussion, we present here a finite size scaling analysis of results from extremely 
long (1/is) MD simulations. We find that 1) LDL is a genuine liquid state, metastable with 
respect to the crystal, 2) LDL and HDL are separated by a first-order phase transition line 
ending in a critical point, 3) LDL has relaxation times that exceed 1/is at temperatures 
below the LDL-HDL coexistence at low pressure, 4) the results are robust with respect to 
the finite size scaling analysis and show that the LDL-HDL critical point belongs to the 3 
dimensional Ising model. 

RESULTS 

To show that the LL phase transition exists in the thermodynamic limit, we perform a 
finite-size analysis along isobars within the supercooled liquid region. For this purpose, we 
calculate the Challa-Landau-Binder par ameter H = 1 — (p^)/3(p^)^ for the bimodality of the 



density distribution function, Vi^p) [20|, l2l|. When V{p) is unimodal, H adopts the value 
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2/3 in the thermodynamic hmit — > oo, while IT < 2/3 when 'P(p) is bimodal, since two 
phases coexist (Fig. [T]d). 

However, for a finite system 11 < 2/3 whenever deviates from a delta function. 

This occurs in the region of the phase diagram where, for a finite system, the isothermal 
compressibility, Kt-, has a maximum, i.e., along a locus in the P-T plane that includes (i) 
the discontinuous (in the thermodynamic limit) phase transition at P > Pc, the LL critical 
pressure, (ii) the effective LL critical point at Pc{N), where the discontinuity vanishes, and 
(iii) a line for P < Pc that emanates from the LL critical point into the supercritical region. 
Near P. this line follows the locus of maxima of the correlation length, known as the Widom 
line 221 , and deviates from it at lower P 23|. 



The finite-size behavior of LI allows us to distinguish whether an isobar is above or below 
\2l\ (Fig. Wp)- When isobars cross the Widom line (P < Pc), LI displays a minimum 



Llniin (inset in Fig. lb) that in leading order approaches 2/3 linearly with 1/A^. When I^(p) 
consists of two Gaussians of equal weight, i.e. at the coexistence line for P > Pc, Llmin 
approaches, also linearly with 1/A^, another limiting value LI — )■ 2 / 3 — ( p"^ — pf)"^ / [3( + pf^] 
where pn = Ph(-P) and pL = Pl(-P) are the densities of the two coexisting phases (20[. This 
limiting value progressively decreases as P increase above Pc, since pu — pl increases at 
coexistence as (P — Pc)'', where /3 ~ 0.3 is the critical exponent of the 3d Ising universality 
class jl?! . 

To ensure that the system is in thermal equilibrium, we calculate the correlation time for 
the first maximum ki of the oxygen-oxygen intermediate scattering function Sooik,t), as 
defined in the Methods Section. While correlation times in the HDL phase are very short 
(^ 0.01 ns), they become of the order of 100 ns in the LDL phase, implying that simulations 
of less than 1 ps are likely affected by poor statistical sampling (Fig. [2]). For temperatures 
above the line Tg in Fig. [1^, correlation times are smaller than 100 ns and we can equilibrate 
the system within our simulation time. 

Figure [3^ shows a typical example of a simulation near the critical point for = 343 
molecules at P = 215 MPa and T = 244 K. Here the system exhibits phase flipping between 
LDL and HDL, with the life-time of each phase distributed from ^ 20 ns to ~ 300ns. This 
nanoscale phase flipping results in a bimodal density distribution (Fig. [Hh") and is observed 
for all temperatures and pressures around the LL phase transition in a region that shrinks 
with growing system size. 
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DISCUSSION 

To estimate the critical exponents of the LL critical point we next investigate the distri- 
bution of the order parameter M of the LL phase transition. As for the liquid-gas phase 



transition 



271], the order parameter is not simply the density, but a linear combination of 



the density with another observable [28|. Here we choose the linear combination of density 
and energy M = p + sE 27j] and find that it follows, as expected, the behavior of a liquid 
in the universality class of the three dimensional (3d) Ising model, as is also the case for the 
liquid-gas transition (Fig. |4]). At P = 205 MPa the difference between the maxima and 
the central minimum of the order parameter distribution is smaller than for the 3d Ising 
case. At P = 210 MPa it is larger and the critical point therefore seems to be in between, 
consistent with the conclusion obtained from the analysis of 11. We get the best fit of the 
order parameter distribution function at a pressure of P = 206 ± 3 MPa and a temperature 
of T = 246 ± 1 K. 

The same analysis for = 512 and 729 yields estimates, consistent with A^ = 343, of the 
LL critical point to be P^ = 208 ± 3 MPa and = 246 ± 1 K (Fig. gb, c). The finite size 
scaling of the amplitudes of the order parameter distribution A ~ L^^'^ is consistent with 
the behavior predicted for the 3d Ising universality class with f3/v ^ 0.518 j^^] and strong 
corrections to scaling for A^ < 343 (Fig. HJi). 

Finally, we investigate also the possibility of spontaneous crystal nucleation in the LDL 
phase using the structural order parameter 19|. At temperatures below the region of 
phase flipping, the samples sometimes form large crystallites filling up to 10% of the system 
volume. Their structure exhibits a mixture of cubic and hexagonal symmetry. However, in 
approximately 99% of simulations these unstable crystallites vanish within the simulation 
time of 1 /is, showing that the free-energy barrier for the crystallization process is signifi- 
cantly larger than ksT in the LDL phase (Fig. |5]). We observe irreversible crystallization in 
only 3 out of 350 (1 /is)-runs, for only A^ = 343 and all corresponding to state points near 
the LL critical point (Fig. [TK). This is consistent with the general result that a metastable 
fluid-fluid phase transition favors the crystallization process in its vicinity 29|. We did not 
observe any crystallization events for A^ = 512 and A^ = 729 although the total simulation 
time for these systems is comparable to that of A^ = 343. The fact that the crystallization 
rate is not increasing with system size is evidence that LDL is the genuine metastable phase 
with respect to the stable crystal phase. 
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In conclusion, we use new methods to investigate both the statics and dynamics of deeply 
supercooled ST2-water. Specifically, we analyze static quantities (density and potential en- 
ergy) using the framework of finite-size scaling theory, and we analyze the dynamic structure 
factor over three orders of magnitude of time scales, from 1 to 1000 ns. We find definitive 
evidence of a first order LL phase transition line between two genuine phases that are each 
metastable with respect to a liquid. The phase transition line terminates in a LL criti- 
cal point, and the exponents associated with this LL critical point are indistinguishable 
from those expected for a three-dimensional lattice-gas model which is used to describe the 
liquid- vapor critical point. 



METHODS 

We performed MD simulations in the NPT ensemble using the Stillinger and Rahman 



34l | five-point water model ST2, consisting of five particles interacting through electrostatic 



and Lennard- Jones forces with a cutoff of 7.8 A. The pressure was not adjusted to correct for 
the effects of the Lennard- Jones cutoff, since it would originate from mean field calculations, 
which become rather poor near a critical point. 

We apply the Shake algorithm to constrain the particles inside each molecule. The con- 
stant pressure is imposed by a Berendsen barostat, and a Nose- Hoover thermostat is applied 
to ensure constant temperature 35|. Periodic boundary conditions have been implemented. 

For the simulations we used the following protocol consisting of three steps: (1) For any 
given density, a constant volume simulation is performed at T = 300 K during 1 ns (first 
pre-run). (2) The ensemble is then changed to NPT by adding the Berendsen barostat with 
the desired pressure and the temperature is reduced to T = 265 K, ensuring that the system 
reaches the HDL phase after 1 ns of equilibration (second pre-run). (3) After these two pre- 
runs the system is quenched to the desired temperature, from which the first 100 — 200 ns 
are removed as thermalization time. The choice of the thermalization time will be discussed 
next. 

To decide whether the equilibration time is sufficient, we perform two steps. First, we 
inspect the time series of energy and density to discard the possibility of spontaneous crys- 
tallization. In all our NPT simulations we observed only three crystallization events (~ 1% 
of total number of runs) all of them in systems with the smallest size [N = 343 molecules). 
We use them as a reference for the crystal. In a second step we measure the correlation time 
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using the intermediate scattering function. 

MD simulations are performed for a finite numbers of state points (Fig. [T^). We use then 
histogram reweighting method jl| to complement the statistics of each state point with the 



information from nearby state points. Histogram reweighting j30| is a method that combines 
the overlapping histograms of quantities calculated at close-enough state points, reweighting 
them with an appropriate factor that takes into account the difference in thermodynamic 
parameters. It is a powerful method that allows to calculate the observables for a continuous 
range of thermodynamic parameters within those directly simulated. 

The order parameter M = p + sE is obtained from the distribution in the density-energy 
plane (Fig. |3^), by integrating it with a delta- function 6{M — p — sE). We select the value 
of s for which the distribution of M best fits the distribution of the order parameter for 
the 3d Ising universality class. The main effect found when changing s is a small shift in 
the estimated critical temperature Tc of about 0.1 K, which is less than the error of 0.5 K 
originating from the histogram reweighting. 

The oxygen-oxygen intermediate scattering function 6*00 (k,t) can be used to distinguish 
between phases of different structure, such as LDL and HDL. We also use it to estimate the 
correlation time. It is defined as 

5oo(k,t) = ^/j]exp(zk-[r,(0-r^(t' + t)])\ , (1) 

\ C,m I ^1 

where denotes averaging over the simulation time t', r^(t') is the position of the oxygen 
of molecule / at time t', k is the wave vector and k is its magnitude |k|. S'oo(k, t) describes 
the time evolution of the spatial correlation along the wave vector k. Since the system has 
periodic boundaries, the components of k have discrete values 27rj/L, where L is the length 
of the simulation box and j = 1,2,.... We define Soo{k,t) = {S{k,t))j, where average is 
taken over all vectors k with magnitude k belonging to jth spherical bin 7r(j — 1/2) /L < 
k < it{j + l/2)/L, for j = 2, 3, ...300. 

The temporal decay of 6*00 (^) t) is characterized by two relaxation times: (i) a short time, 
Tg, after which Soo{k, t) reaches a plateau Soo{k, Tp) corresponding to the bouncing of the 
particles inside the cages formed by their neighbors, and (ii) a long time, Tq,, corresponding 
to a particle escaping from its cage and diffusing away from its initial position. We define 
the correlation time the time for which Coo{k,T) = Soo{k,T)/ Soo{k,Ti3) = 1/e, 

where Coo{k, r) is the structural correlation function (Fig. [2]). 



We define the bond order parameter ^3 following Ref. 19|. The quantity d-i{i,j) charac- 
terizes the bond between molecules i and j and is designed to distinguish between a fluid 
and a diamond structure. It uses the Y-^ spherical harmonics to identify the tetragonal 
symmetry of the diamond structure. In general, each molecule is characterized by a vector 

in the (4£ + 2)-dimensional Euclidean space with components TZe{q\^ and Xm{q\^ 
{m = -1,0, !,....,£), with 



If molecule j belongs to the first coordination shell rij (shell of four nearest neighbors) of 
molecule we define ds{i,j) as the cosine of the angle between two vectors and q3 
characterizing the first coordination shells of molecules j and i, respectively: 

m\m\ 

where 

(qa • qD = Yl li^'^^ iim + ^^ 4,m^^ qIJ, 

m=-l 

and |q*3| = V(q^qi). 

In a perfect diamond crystal d3{i,j) = —1 for all bonds, while for a graphite crystal 
dzihi) = ^1 only for bonds connecting atoms in the same layer. For bonds connecting 
atoms in different layers d^{i^j) = —1/9. Thus in graphite each atom has three out of four 
bonds having ds{i,j) = —1. In our simulations, the spontaneously grown crystals have many 
defects, with different parts of the crystals following diamond or graphite patterns (Fig. [6|l). 
Therefore, we consider a molecule in a crystal to have either three or four bonds with 
dzihi) < dc = —0.87, where the value of dc = —0.87 is selected as two standard deviations 
from the peak of the crystal histogram corresponding to d^ = —1. This is exactly the same 
criterion to specify molecules in the crystal state as in Ref. 19[ . We find separate crystallites 
using the percolation criterion, i.e., two molecules satisfying the crystalline criterion belong 
to the same crystallite if they belong to the first coordination shell of each other (Fig. |5]). 

We finally observe that by defining ip3{i) = j J2'j=i d^ihi) as the average of ^3 over the 
four bonds of each molecule, we introduce a single-molecule structural parameter that also 
can be used to distinguish among the HDL, the LDL and the crystal phase (Fig. Wp)- 
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FIGURE LEGENDS 

FIG. 1: Phase diagram and finite size scaling analysis to locate the line of liquid- 
liquid (LL) phase transitions, (a) State points in the P-T diagram simulated. Different 
symbols correspond to different sizes N. The high-T (red) region exhibits HDL-like states and the 
low-T (blue) region LDL-like states. In the intermediate (violet) region we observe flipping between 
HDL-like and LDL-like states. Below the black line correlation times are larger than 100 ns, while 
above they are smaller. Equilibrium is attained within reasonable simulation times. The white 
region, denoted CP, is our estimate of the location of the LL critical point in the thermodynamic 
limit, (b) Finite-size analysis of Ilniin along isobars crossing the discontinuous LL phase transition 
(violet at high P) and the Widom line (within the violet region at low P). At P = 190 MPa, Iljnin 
approaches 2/3 when N ^ oo, indicating that the density distribution is unimodal and that one 
crosses the Widom line, and not the line of discontinuous phase transition. At P = 200 MPa, Ilmin 
app roaches ~ 2/3 — 0.001, consistent within its error bar with the value expected at coexistence 
20|. At P = 210 MPa, Hmm tends to a smaller value clearly excluding 2/3 and therefore the 
distribution 'D{p) is bimodal, that is the fingerprint of a discontinuous LL phase transition. Ilmin 
depends linearly on 1/A^ to the leading order, displaying deviations only for the smallest size 
N = 216. The inset shows 11 along the isobar at P = 200 MPa as a function of T for all four 
system sizes (from bottom to top: N = 216, 343, 512, 729) displaying a clear minimum Ilmin. Lines 
are interpolations obtained using histogram reweighting for up to eleven independent simulations 
of length 1 /i. 
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FIG. 2: Definition of the correlation time tq using the intermediate scattering function. 

The correlation time tq is calculated using the correlation function Cooik,t) of the intermediate 
scattering function of the oxygen atoms Soo{k,t). For the k vectors corresponding to the first 
three maxima ki, k2 and k^ (marked in red, blue and green in the inset), we calculate the evolution 
of the correlation function Coo{ki,T). We then define the correlation time as the time for which 
Cooiki,T) decreases to 1/e for the slowest of the ki vectors. For nearly all the state points ki has 
been the vector for which this decrease has been the slowest. 10-100 ns for the LDL phase, so we 
can equilibrate this phase in our simulations of about 1000 ns. Data are for a system of = 343 
molecules at pressure P = 210 MPa and temperatures (from left to right) T = 244 K, 243 K, in 
the LDL phase, and 242 K below the Tg line of Fig. [1^. 

FIG. 3: Phase flipping between LDL and HDL at coexistence, (a) The 1 fis time series 
shows how frequently, at constant P = 215 MPa and T = 244 K, = 343 ST2- water molecules 
switch from LDL-like to HDL-like states, (b) The histogram for the sampled density values, in 
arbitrary units, after discarding the first 100 ns of the 1 /is time series. For LDL-like states 
p ~ (0.89 lb 0.01) g/cm^ and for HDL-like states p ~ (1.02 it 0.03) g/cm^ corresponding to a 
difference of ~ 13% in density. Dashed lines are Gaussian best fits of the histogram around the 
two maxima. 
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FIG. 4: The liquid-liquid critical point falls into the same universality class as the liquid- 
gas critical point, (a) The distribution function of the rescaled order parameter x = A{M — Mc) 
where M = p + sE with s = 27.60^^, follows for P = (206 ± 3) MPa and T = (246 ± 1) K 
(triangles) the order parameter distribution function of the 3d Ising model (solid line) [3a] . The data 
are from histogram reweighting of = 343 molecules at P = 205 MPa and T = 246.6 K (squares), 
P = 206 MPa and T = 246 K (triangles) and P = 210 MPa and T = 245.1 K (circles). We repeat 
the analysis for (b) N = 512 and (c) N = 729. (d) For large sizes the amplitude A (triangles) scales 
as ^ L^/*^, where /3/i^ ~ 0.52, as in the 3d Ising universality class For N < 343 corrections 
to scaling are strong, (e) Contour plot of the distribution of states in the density-energy plane, 
with red corresponding to the highest values and blue to the lowest. The distribution of the order 
parameter M = p + sE is obtained from this two-dimensional distribution by integrating it with 
a delta-function 6{M — p — sE). We select the value of s for which the distribution of M best fits 
the distribution of the order parameter for the 3d Ising universality class. 



FIG. 5: Example of a simulation where the largest crystallite grows up to 35 molecules 
and then vanishes in a system having N = 343 molecules at a pressure of P = 200 MPa 

and T = 246 K. A molecule i is considered to belong to a crystal if d^{i,j) < dc = —0.87 for three 
out of its four bonds with nearest neighbors j. 
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FIG. 6: The equilibrium probability distributions of two different structural parameters 
that allow us to distinguish among three structures: HDL, LDL, and crystal, (a) For 

the smallest size simulated {N = 343), in the vicinity of the LL critical point in ~ 1% of the runs 
our system spontaneously crystallizes, forming a structures with diamond and graphite patterns 
and with many defects. This structure is well characterized by the probability distribution V^ds) 
(line with full dots) of the parameter ^3, displaying a large maxima close to ^3 = —1, the value that 
corresponds to the perfect diamond crystal. For the sake of comparison with the other cases, we 
divide V^d^) of the crystal by 10. In the 99% of our simulations we find distributions '^'(^3) as those 
presented here for P = 215 MPa, a pressure above the LL crital point pressure, and decreasing T 
(from right to left). ^(^3) shows an abrupt change when crossing the first-order LL phase transition 
region at T ~ 244 K. In particular, 'D{d3) displays a pronounced shoulder at higher ^3 in the LDL 
phase, and is very different from the crystal case. The arrow marks the value dc = —0.87 selected 
as two standard deviations from the peak of the crystal histogram corresponding to ^3 = — 1 as 
in Ref. [l^. (b) The equilibrium probability distribution of the single-molecule parameter ip^ also 
distinguishes among HDL, LDL and crystal structures. The average values for fluid phases are 
= -0.34 ± 0.19 in the HDL phase and Tps = -0.57 ± 0.16 in the LDL phase. 
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